Honeybees are very important organisms to the world at large. Honeybees are important pollinators of flowering plants, and thus are essential to the agriculture industry. Despite their vital importance, honeybee populations have been declining for nearly thirty years. One of the main reasons for this decline in honeybee numbers is the Varroa mite (Varroa destructor), an external parasite of honeybees. Varroa mites are responsible for causing Varroosis in honeybees, weakening adult bees by feeding on their haemolymph and thereby increasing their susceptibility to other diseases. The current protocol in place to help rid honeybee colonies of Varroa mites is through miticides. However, this method of pest control is not sustainable as the use of miticides may further weaken the honeybees or contaminate hive products. More importantly, an increase in resistance to miticides by Varroa mites has been seen. Consequently, it is imperative to develop a new means of controlling these pests to increase honeybee populations in the long term.
This research aims to examine various odorants and their effects on both Varroa mites and honeybees to analyze their attractant or repellent properties on the two species. Additionally, the effects of the odorants on honeybee and Varroa mite behaviour, as well as effects on in-hive production, will be examined.
The outcome of this research project will be to implement the use of odorants as a more sustainable means of controlling Varroa mite populations, without the use of toxic pesticides. Successful implementation of the results of this study could help slow the decline of honeybee populations worldwide by decreasing the presence of Varroa mites in honeybee colonies.
Over the summer, Mike conducted electrotarsograms (ETGs) on Varroa mites to examine their physiological sensitivity to several previously identified attractive odourants.
ETGs on Varroa mites were conducted using single odourants previously identified as evoking a response in Varroa mites, and were tested in four different concentrations (10ng, 1ng, 100μg, 10μg). These are typical concentrations used in ETG studies and often encompass concentrations close to what Varroa mites would encounter in their natural environment. Relative to a control stimulus, these reactions can be used to infer the relative importance of each odourant to Varroa mites in the hive environment.
The questions we want to answer from the dataset are as follows:
1. Are there mite trials that are outliers?
2. Are there odours which consistantly evoke a response?
3. Are there concentration-dependent trends?
4. Can we quantify and account for between-trial (between-mite) variability?
my.df <- read_csv('DataForR.csv')
## Parsed with column specification:
## cols(
## trial = col_integer(),
## odour = col_character(),
## conc = col_double(),
## percent = col_double(),
## gb = col_character(),
## normalized = col_double()
## )
str(my.df)
## Classes 'tbl_df', 'tbl' and 'data.frame': 1165 obs. of 6 variables:
## $ trial : int 1 1 1 1 1 1 1 1 1 1 ...
## $ odour : chr "2-Hydroxyhexanoic Acid (DCM)" "2-Hydroxyhexanoic Acid (DCM)" "2-Hydroxyhexanoic Acid (DCM)" "Methyl Palmitate" ...
## $ conc : num 0.01 0.1 1 0.01 0 10 0.01 1 0.01 0.1 ...
## $ percent : num 53.6 142.5 157.9 137.7 100 ...
## $ gb : chr "Good" "Good" "Good" "Good" ...
## $ normalized: num 13.3 129.8 136.7 127.4 100 ...
## - attr(*, "spec")=List of 2
## ..$ cols :List of 6
## .. ..$ trial : list()
## .. .. ..- attr(*, "class")= chr "collector_integer" "collector"
## .. ..$ odour : list()
## .. .. ..- attr(*, "class")= chr "collector_character" "collector"
## .. ..$ conc : list()
## .. .. ..- attr(*, "class")= chr "collector_double" "collector"
## .. ..$ percent : list()
## .. .. ..- attr(*, "class")= chr "collector_double" "collector"
## .. ..$ gb : list()
## .. .. ..- attr(*, "class")= chr "collector_character" "collector"
## .. ..$ normalized: list()
## .. .. ..- attr(*, "class")= chr "collector_double" "collector"
## ..$ default: list()
## .. ..- attr(*, "class")= chr "collector_guess" "collector"
## ..- attr(*, "class")= chr "col_spec"
summary(my.df)
## trial odour conc percent
## Min. : 1.00 Length:1165 Min. : 0.000 Min. : 17.59
## 1st Qu.: 9.00 Class :character 1st Qu.: 0.010 1st Qu.: 78.07
## Median :13.00 Mode :character Median : 0.100 Median : 100.00
## Mean :12.22 Mean : 2.665 Mean : 117.17
## 3rd Qu.:16.00 3rd Qu.: 1.000 3rd Qu.: 125.02
## Max. :18.00 Max. :10.000 Max. :2507.58
## gb normalized
## Length:1165 Min. :-368.60
## Class :character 1st Qu.: 71.91
## Mode :character Median : 100.00
## Mean : 89.80
## 3rd Qu.: 120.01
## Max. : 196.01
hist(my.df$normalized,
xlab="Normalized ETG Response", ylab="Frequency",
main="Frequency of Normalized ETG Responses")
p <- ggplot(data=my.df, aes(normalized))
p + geom_density()
Doing this produces a histogram of the percent reaction by Varroa mites. This is an easy way to see how the data is distributed.
#change integer to factor for trials
my.df$trial<-factor(my.df$trial)
#lists the levels of factors
levels(my.df$trial)
## [1] "1" "2" "4" "5" "6" "7" "8" "9" "10" "11" "12" "13" "14" "15"
## [15] "16" "17" "18"
#change number to factor for concentration "conc"
my.df$conc<-factor(my.df$conc)
levels(my.df$conc)
## [1] "0" "0.01" "0.1" "1" "10"
# types of odours
my.df$odour<-factor(my.df$odour)
levels(my.df$odour)
## [1] "1-Hexadecanol" "2-Heptanol"
## [3] "2-Heptanone" "2-Hydroxyhexanoic Acid (DCM)"
## [5] "2-Nonanal" "Benzoic Acid"
## [7] "Butyric Acid" "Dodecyl Aldehyde"
## [9] "Ethyl Palmitate" "Geraniol"
## [11] "Heptacosane" "Heptadecane"
## [13] "Hexadecanal" "Hexane"
## [15] "Methyl Oleate" "Methyl Palmitate"
## [17] "Octadecanol" "Octanoic Acid"
## [19] "Palmitic Acid" "trans-Nerolidol"
We discovered that these runs were highly variable and produced messy plots, so we will remove them for now.
good_runs<- filter(my.df,
!(odour %in% c("1-Hexadecanol", "Methyl Palmitate", "Hexane", "Palmitic Acid" )))
subsetted <- mutate(good_runs,
Zresponse = normalized - 100)
hist(subsetted$Zresponse,xlab="Normalized ETG Response",ylab="Frequency", main="Frequency of Normalized ETG Responses")
A negative response makes no sense physiologically, so this is an indication of a bad run. We remove these negative values from analysis.
nonegatives <- subsetted %>%
mutate(Zresponse= ifelse(Zresponse <0, NA, Zresponse))
str(nonegatives)
## Classes 'tbl_df', 'tbl' and 'data.frame': 706 obs. of 7 variables:
## $ trial : Factor w/ 17 levels "1","2","4","5",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ odour : Factor w/ 20 levels "1-Hexadecanol",..: 4 4 4 15 12 12 12 2 5 5 ...
## $ conc : Factor w/ 5 levels "0","0.01","0.1",..: 2 3 4 5 2 3 4 3 2 3 ...
## $ percent : num 53.6 142.5 157.9 202.7 44.1 ...
## $ gb : chr "Good" "Good" "Good" "Good" ...
## $ normalized: num 13.3 129.8 136.7 150.7 -27 ...
## $ Zresponse : num NA 29.8 36.7 50.7 NA ...
Now that we have manipulated the data in a way that is useful to us, we will make our faceted plot.
ggplot(nonegatives)+
geom_boxplot(aes(x=conc,y= Zresponse, group= conc))+
xlab("Concentration")+ylab("Normalized Percent Reaction")+ggtitle("Faceted ETG Plots")+
theme_bw(10)+
facet_wrap(~odour)
## Warning: Removed 366 rows containing non-finite values (stat_boxplot).
From this plot, we can see the overall trend in ETG responses from Varroa mites at each concentration for each compound tested.
We are still waiting on temperature and humidity data, but Mike has been in touch with someone regarding this.
absoluteweather<- read.csv('absoluteweather.csv')
library(tidyverse)
absoluteweather$conc <- factor(absoluteweather$conc)
absoluteweather$trial <- factor(absoluteweather$trial)
str(absoluteweather)
## 'data.frame': 1165 obs. of 13 variables:
## $ trial : Factor w/ 17 levels "1","2","4","5",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ odour : Factor w/ 20 levels "1-Hexadecanol",..: 4 4 4 16 14 15 1 16 12 12 ...
## $ conc : Factor w/ 5 levels "0","0.01","0.1",..: 2 3 4 2 1 5 2 4 2 3 ...
## $ absamp : num 3.02 8.09 9 7.95 5.91 ...
## $ linearhex : num 5.64 5.67 5.7 5.77 5.91 ...
## $ percent : num 53.6 142.5 157.9 137.7 100 ...
## $ normalized: num 13.3 129.8 136.7 127.4 100 ...
## $ gb : Factor w/ 2 levels "Bad","Good": 2 2 2 2 2 2 2 2 2 2 ...
## $ humid : num 69 69 69.1 69.2 69.3 ...
## $ temp : num 20.8 20.8 20.8 20.7 20.7 ...
## $ dewp : num 14.9 14.9 14.9 14.9 14.9 ...
## $ hour : num 2013 2014 2014 2016 2019 ...
## $ min : num 13.1 13.8 14.5 16 18.9 ...
ggplot(data = absoluteweather, mapping = aes(humid, absamp)) +
geom_jitter(aes(colour = absoluteweather$trial), width = 0.25) +
ggtitle("Amplitude vs Humidity")
ggplot(data = absoluteweather, mapping = aes(temp, absamp)) +
geom_jitter(aes(colour = absoluteweather$trial), width = 0.25) +
ggtitle("Amplitude vs Temperature")
ggplot(data = absoluteweather, mapping = aes(dewp, absamp)) +
geom_jitter(aes(colour = absoluteweather$trial), width = 0.25) +
ggtitle("Amplitude vs Dewpoint")
ggplot(data = absoluteweather, mapping = aes(absamp))+
geom_histogram(bins = 50)
# correlations
cor(absoluteweather$absamp, absoluteweather$temp)
## [1] -0.4868699
cov(absoluteweather$absamp, absoluteweather$temp)
## [1] -234.6016
corr_amp_tempr <- cor.test(x=absoluteweather$absamp,
y=absoluteweather$temp, method = 'spearman')
## Warning in cor.test.default(x = absoluteweather$absamp, y = absoluteweather
## $temp, : Cannot compute exact p-value with ties
corr_amp_tempr
##
## Spearman's rank correlation rho
##
## data: absoluteweather$absamp and absoluteweather$temp
## S = 367540000, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## -0.3946761
corr_amp_humid <- cor.test(x=absoluteweather$absamp,
y=absoluteweather$humid, method = 'spearman')
## Warning in cor.test.default(x = absoluteweather$absamp, y = absoluteweather
## $humid, : Cannot compute exact p-value with ties
corr_amp_humid
##
## Spearman's rank correlation rho
##
## data: absoluteweather$absamp and absoluteweather$humid
## S = 103650000, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.606672
corr_amp_dewp<- cor.test(x=absoluteweather$absamp,
y=absoluteweather$dewp, method = 'spearman')
## Warning in cor.test.default(x = absoluteweather$absamp, y = absoluteweather
## $dewp, : Cannot compute exact p-value with ties
corr_amp_dewp
##
## Spearman's rank correlation rho
##
## data: absoluteweather$absamp and absoluteweather$dewp
## S = 195510000, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.2580989
library(fifer)
## Loading required package: MASS
##
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
##
## select
spearman.plot(absoluteweather$absamp, absoluteweather$temp)
weather_matrix <- data.matrix(absoluteweather[7:10], rownames.force = NA)
pairs(weather_matrix) # plots all the stuff simultaneiously
absoluteweather <- read.csv('absoluteweather.csv')
str(absoluteweather)
## 'data.frame': 1165 obs. of 13 variables:
## $ trial : int 1 1 1 1 1 1 1 1 1 1 ...
## $ odour : Factor w/ 20 levels "1-Hexadecanol",..: 4 4 4 16 14 15 1 16 12 12 ...
## $ conc : num 0.01 0.1 1 0.01 0 10 0.01 1 0.01 0.1 ...
## $ absamp : num 3.02 8.09 9 7.95 5.91 ...
## $ linearhex : num 5.64 5.67 5.7 5.77 5.91 ...
## $ percent : num 53.6 142.5 157.9 137.7 100 ...
## $ normalized: num 13.3 129.8 136.7 127.4 100 ...
## $ gb : Factor w/ 2 levels "Bad","Good": 2 2 2 2 2 2 2 2 2 2 ...
## $ humid : num 69 69 69.1 69.2 69.3 ...
## $ temp : num 20.8 20.8 20.8 20.7 20.7 ...
## $ dewp : num 14.9 14.9 14.9 14.9 14.9 ...
## $ hour : num 2013 2014 2014 2016 2019 ...
## $ min : num 13.1 13.8 14.5 16 18.9 ...
absoluteweather$trial <- as.factor(absoluteweather$trial)
absoluteweather$conc <- as.factor(absoluteweather$conc)
subsetted <- mutate(absoluteweather,
Zresponse = normalized - 100)
#Not sure whether to replace nonsense scores with NA or 0
nonegatives <- subsetted %>%
mutate(Zresponse= ifelse(Zresponse <0, NA, Zresponse))
str(nonegatives$Zresponse)
## num [1:1165] NA 29.8 36.7 27.4 0 ...
ggplot(data = nonegatives, aes(min, absamp))+
geom_jitter(aes(colour = nonegatives$odour), width = 0.25)+
xlab("Time")+ylab("Absamp")+ggtitle("Faceted ETG Plots")+
theme_bw(10)+
facet_wrap(~trial)
nonegatives$trial <- as.numeric(nonegatives$trial)
goodtrials <- nonegatives[nonegatives$trial %in% c(1:8,11, 12), ]
goodtrials$trial <- as.factor(goodtrials$trial)
str(goodtrials)
## 'data.frame': 530 obs. of 14 variables:
## $ trial : Factor w/ 10 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ odour : Factor w/ 20 levels "1-Hexadecanol",..: 4 4 4 16 14 15 1 16 12 12 ...
## $ conc : Factor w/ 5 levels "0","0.01","0.1",..: 2 3 4 2 1 5 2 4 2 3 ...
## $ absamp : num 3.02 8.09 9 7.95 5.91 ...
## $ linearhex : num 5.64 5.67 5.7 5.77 5.91 ...
## $ percent : num 53.6 142.5 157.9 137.7 100 ...
## $ normalized: num 13.3 129.8 136.7 127.4 100 ...
## $ gb : Factor w/ 2 levels "Bad","Good": 2 2 2 2 2 2 2 2 2 2 ...
## $ humid : num 69 69 69.1 69.2 69.3 ...
## $ temp : num 20.8 20.8 20.8 20.7 20.7 ...
## $ dewp : num 14.9 14.9 14.9 14.9 14.9 ...
## $ hour : num 2013 2014 2014 2016 2019 ...
## $ min : num 13.1 13.8 14.5 16 18.9 ...
## $ Zresponse : num NA 29.8 36.7 27.4 0 ...
remove(subsetted)
ggplot(data = nonegatives, mapping = aes(Zresponse))+
geom_histogram(bins = 50)
## Warning: Removed 577 rows containing non-finite values (stat_bin).
plot(aov_out <- aov(Zresponse~odour*conc*trial, data = goodtrials))
## Warning: not plotting observations with leverage one:
## 1, 2, 5, 8, 12, 13, 14, 19, 21, 22, 23, 24, 26, 28, 31, 32, 33, 38, 39, 40, 43, 45, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 64, 65, 67, 68, 69, 70, 71, 85, 86, 87, 88, 89, 90, 96, 97, 98, 99, 100, 101, 102, 103, 104, 105, 106, 107, 116, 117, 118, 119, 120, 125, 126, 127, 128, 129, 131, 132, 133, 136, 137, 140, 141, 146, 154, 157, 158, 159, 160, 161, 162, 163, 164, 165, 166, 171, 172, 177, 178, 179, 180, 181, 183, 184, 185, 186, 187, 191, 192, 193, 194, 195, 196, 197, 198, 199, 205, 206, 207, 208, 209, 210, 211, 212, 217, 218, 219, 222, 230, 233, 238, 239, 240, 241, 242, 243, 244, 251, 252, 253, 254, 257, 258, 260, 267, 268, 269
## Warning: not plotting observations with leverage one:
## 1, 2, 5, 8, 12, 13, 14, 19, 21, 22, 23, 24, 26, 28, 31, 32, 33, 38, 39, 40, 43, 45, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 64, 65, 67, 68, 69, 70, 71, 85, 86, 87, 88, 89, 90, 96, 97, 98, 99, 100, 101, 102, 103, 104, 105, 106, 107, 116, 117, 118, 119, 120, 125, 126, 127, 128, 129, 131, 132, 133, 136, 137, 140, 141, 146, 154, 157, 158, 159, 160, 161, 162, 163, 164, 165, 166, 171, 172, 177, 178, 179, 180, 181, 183, 184, 185, 186, 187, 191, 192, 193, 194, 195, 196, 197, 198, 199, 205, 206, 207, 208, 209, 210, 211, 212, 217, 218, 219, 222, 230, 233, 238, 239, 240, 241, 242, 243, 244, 251, 252, 253, 254, 257, 258, 260, 267, 268, 269
## Warning in sqrt(crit * p * (1 - hh)/hh): NaNs produced
## Warning in sqrt(crit * p * (1 - hh)/hh): NaNs produced
summary(aov_out)
## Df Sum Sq Mean Sq F value Pr(>F)
## odour 19 39802 2095 6.955 8.29e-11 ***
## conc 3 523 174 0.579 0.6305
## trial 9 40810 4534 15.054 1.64e-14 ***
## odour:conc 41 11033 269 0.893 0.6497
## odour:trial 55 23456 426 1.416 0.0722 .
## conc:trial 25 6152 246 0.817 0.7108
## odour:conc:trial 31 7147 231 0.765 0.7978
## Residuals 88 26507 301
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 258 observations deleted due to missingness
library(tidyverse)
library(modelr)
options(na.action = na.warn)
ggplot(goodtrials, aes(trial, absamp))+
geom_point(aes(colour=odour))
Zresponse ~ concentration+odour
mod01 <- lm(Zresponse~conc+odour, data = goodtrials)
## Warning: Dropping 258 rows with missing values
#fits red dots as predictions according to the model
grid <- goodtrials %>%
data_grid(odour,conc) %>%
add_predictions(mod01)
## Warning in predict.lm(model, data): prediction from a rank-deficient fit
## may be misleading
grid
## # A tibble: 100 x 3
## odour conc pred
## <fct> <fct> <dbl>
## 1 1-Hexadecanol 0 0.000000000000165
## 2 1-Hexadecanol 0.01 27.2
## 3 1-Hexadecanol 0.1 29.9
## 4 1-Hexadecanol 1 30.9
## 5 1-Hexadecanol 10 27.6
## 6 2-Heptanol 0 - 7.80
## 7 2-Heptanol 0.01 19.4
## 8 2-Heptanol 0.1 22.1
## 9 2-Heptanol 1 23.1
## 10 2-Heptanol 10 19.8
## # ... with 90 more rows
#plots our data
ggplot(goodtrials, aes(x = conc)) +
geom_point(aes(y = Zresponse)) +
geom_point(data = grid, aes(y = pred), colour = "red", size = 0.2)+
facet_wrap (~odour)
## Warning: Removed 258 rows containing missing values (geom_point).
This model doesn’t appear to fit the data very well
#I tells the below function to change Zresponse to an integer
mod01 <- glm(I(as.integer(goodtrials$Zresponse))~conc+odour, data = goodtrials, poisson)
## Warning: Dropping 258 rows with missing values
## Warning: glm.fit: algorithm did not converge
par(mfrow=c(2, 2))
summary(mod01)
##
## Call:
## glm(formula = I(as.integer(goodtrials$Zresponse)) ~ conc + odour,
## family = poisson, data = goodtrials)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -41.048 -5.423 -1.785 1.715 9.571
##
## Coefficients: (1 not defined because of singularities)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 6.736328 0.003304 2038.804 < 2e-16
## conc0.01 -3.454085 0.095869 -36.029 < 2e-16
## conc0.1 -3.363648 0.100125 -33.594 < 2e-16
## conc1 -3.332014 0.098638 -33.780 < 2e-16
## conc10 -3.442636 0.100819 -34.147 < 2e-16
## odour2-Heptanol -0.307637 0.118917 -2.587 0.009682
## odour2-Heptanone 0.234648 0.107940 2.174 0.029715
## odour2-Hydroxyhexanoic Acid (DCM) -0.080565 0.111600 -0.722 0.470348
## odour2-Nonanal 0.024033 0.109300 0.220 0.825966
## odourBenzoic Acid -0.608394 0.123951 -4.908 9.18e-07
## odourButyric Acid -0.642706 0.178872 -3.593 0.000327
## odourDodecyl Aldehyde 0.050435 0.127273 0.396 0.691904
## odourEthyl Palmitate 0.265458 0.138084 1.922 0.054550
## odourGeraniol -0.785819 0.148409 -5.295 1.19e-07
## odourHeptacosane -0.130434 0.126651 -1.030 0.303071
## odourHeptadecane 0.185362 0.106257 1.744 0.081079
## odourHexadecanal 0.278775 0.110403 2.525 0.011568
## odourHexane NA NA NA NA
## odourMethyl Oleate 0.274981 0.106797 2.575 0.010030
## odourMethyl Palmitate 0.116862 0.098339 1.188 0.234691
## odourOctadecanol -0.198848 0.157146 -1.265 0.205739
## odourOctanoic Acid 0.034727 0.110097 0.315 0.752441
## odourPalmitic Acid 0.598129 0.174819 3.421 0.000623
## odourtrans-Nerolidol -0.231625 0.177177 -1.307 0.191108
##
## (Intercept) ***
## conc0.01 ***
## conc0.1 ***
## conc1 ***
## conc10 ***
## odour2-Heptanol **
## odour2-Heptanone *
## odour2-Hydroxyhexanoic Acid (DCM)
## odour2-Nonanal
## odourBenzoic Acid ***
## odourButyric Acid ***
## odourDodecyl Aldehyde
## odourEthyl Palmitate .
## odourGeraniol ***
## odourHeptacosane
## odourHeptadecane .
## odourHexadecanal *
## odourHexane
## odourMethyl Oleate *
## odourMethyl Palmitate
## odourOctadecanol
## odourOctanoic Acid
## odourPalmitic Acid ***
## odourtrans-Nerolidol
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 6453.5 on 271 degrees of freedom
## Residual deviance: 71316.0 on 249 degrees of freedom
## (258 observations deleted due to missingness)
## AIC: 72458
##
## Number of Fisher Scoring iterations: 25
plot(mod01)
## Warning: not plotting observations with leverage one:
## 520
## Warning: not plotting observations with leverage one:
## 520
This model doesn’t appear to fit the data well
mod02 <- glm(goodtrials$Zresponse~conc+odour, data = goodtrials)
## Warning: Dropping 258 rows with missing values
par(mfrow=c(2, 2))
summary(mod02)
##
## Call:
## glm(formula = goodtrials$Zresponse ~ conc + odour, data = goodtrials)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -36.056 -14.001 -0.068 9.043 65.476
##
## Coefficients: (1 not defined because of singularities)
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.654e-13 3.400e+00 0.000 1.00000
## conc0.01 2.724e+01 1.132e+01 2.406 0.01684
## conc0.1 2.993e+01 1.184e+01 2.528 0.01209
## conc1 3.091e+01 1.168e+01 2.647 0.00864
## conc10 2.757e+01 1.189e+01 2.320 0.02118
## odour2-Heptanol -7.799e+00 1.290e+01 -0.604 0.54613
## odour2-Heptanone 7.476e+00 1.254e+01 0.596 0.55152
## odour2-Hydroxyhexanoic Acid (DCM) -2.186e+00 1.250e+01 -0.175 0.86132
## odour2-Nonanal 5.983e-01 1.240e+01 0.048 0.96154
## odourBenzoic Acid -1.290e+01 1.273e+01 -1.014 0.31172
## odourButyric Acid -1.326e+01 1.650e+01 -0.804 0.42219
## odourDodecyl Aldehyde 1.108e+00 1.457e+01 0.076 0.93942
## odourEthyl Palmitate 7.835e+00 1.660e+01 0.472 0.63739
## odourGeraniol -1.615e+01 1.410e+01 -1.145 0.25323
## odourHeptacosane -3.669e+00 1.397e+01 -0.263 0.79311
## odourHeptadecane 5.537e+00 1.223e+01 0.453 0.65110
## odourHexadecanal 8.822e+00 1.297e+01 0.680 0.49700
## odourHexane NA NA NA NA
## odourMethyl Oleate 8.774e+00 1.243e+01 0.706 0.48102
## odourMethyl Palmitate 3.294e+00 1.113e+01 0.296 0.76752
## odourOctadecanol -5.052e+00 1.660e+01 -0.304 0.76114
## odourOctanoic Acid 7.602e-01 1.254e+01 0.061 0.95169
## odourPalmitic Acid 2.171e+01 2.433e+01 0.892 0.37322
## odourtrans-Nerolidol -6.257e+00 1.865e+01 -0.336 0.73750
##
## (Intercept)
## conc0.01 *
## conc0.1 *
## conc1 **
## conc10 *
## odour2-Heptanol
## odour2-Heptanone
## odour2-Hydroxyhexanoic Acid (DCM)
## odour2-Nonanal
## odourBenzoic Acid
## odourButyric Acid
## odourDodecyl Aldehyde
## odourEthyl Palmitate
## odourGeraniol
## odourHeptacosane
## odourHeptadecane
## odourHexadecanal
## odourHexane
## odourMethyl Oleate
## odourMethyl Palmitate
## odourOctadecanol
## odourOctanoic Acid
## odourPalmitic Acid
## odourtrans-Nerolidol
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 462.2682)
##
## Null deviance: 155429 on 271 degrees of freedom
## Residual deviance: 115105 on 249 degrees of freedom
## (258 observations deleted due to missingness)
## AIC: 2464.9
##
## Number of Fisher Scoring iterations: 2
plot(mod02)
## Warning: not plotting observations with leverage one:
## 520
## Warning: not plotting observations with leverage one:
## 520
This model appears to fit the data better than the poisson model
bindf <- goodtrials %>%
mutate(BinZresponse= ifelse(Zresponse > 0, 1, 0))
mod03 <- glm(bindf$BinZresponse~conc+odour, data = bindf, binomial)
## Warning: Dropping 258 rows with missing values
## Warning: glm.fit: algorithm did not converge
par(mfrow=c(2, 2))
summary(mod03)
##
## Call:
## glm(formula = bindf$BinZresponse ~ conc + odour, family = binomial,
## data = bindf)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.409e-06 2.409e-06 2.409e-06 2.409e-06 2.409e-06
##
## Coefficients: (1 not defined because of singularities)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.657e+01 5.631e+04 0 1
## conc0.01 5.313e+01 1.875e+05 0 1
## conc0.1 5.313e+01 1.961e+05 0 1
## conc1 5.313e+01 1.935e+05 0 1
## conc10 5.313e+01 1.969e+05 0 1
## odour2-Heptanol 2.117e-08 2.137e+05 0 1
## odour2-Heptanone 2.471e-09 2.077e+05 0 1
## odour2-Hydroxyhexanoic Acid (DCM) 6.009e-09 2.071e+05 0 1
## odour2-Nonanal 2.924e-09 2.053e+05 0 1
## odourBenzoic Acid -2.023e-09 2.109e+05 0 1
## odourButyric Acid -4.927e-06 2.732e+05 0 1
## odourDodecyl Aldehyde 9.531e-10 2.413e+05 0 1
## odourEthyl Palmitate -3.604e-08 2.750e+05 0 1
## odourGeraniol 3.838e-08 2.336e+05 0 1
## odourHeptacosane -9.987e-09 2.314e+05 0 1
## odourHeptadecane -1.393e-09 2.025e+05 0 1
## odourHexadecanal -4.719e-10 2.148e+05 0 1
## odourHexane NA NA NA NA
## odourMethyl Oleate 1.646e-09 2.059e+05 0 1
## odourMethyl Palmitate 2.710e-09 1.844e+05 0 1
## odourOctadecanol -3.604e-08 2.750e+05 0 1
## odourOctanoic Acid -4.931e-06 2.076e+05 0 1
## odourPalmitic Acid -1.437e-08 4.030e+05 0 1
## odourtrans-Nerolidol 2.765e-08 3.089e+05 0 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 2.2716e+02 on 271 degrees of freedom
## Residual deviance: 1.5780e-09 on 249 degrees of freedom
## (258 observations deleted due to missingness)
## AIC: 46
##
## Number of Fisher Scoring iterations: 25
plot(mod03)
## Warning: not plotting observations with leverage one:
## 520
## Warning: not plotting observations with leverage one:
## 520